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Abstract 



Within systems biology there is an increasing interest in the stochastic behavior of genetic 
and biochemical reaction networks. An appropriate stochastic description is provided by the 
chemical master equation, which represents a continuous time Markov chain (CTMC). In 
this paper we consider the stochastic properties of a biochemical circuit, known to control 
eukaryotic cell cycle and possibly involved in oncogenesis, recently proposed in the literature 
within a deterministic framework. Due to the inherent stochasticity of biochemical processes 
and the small number of molecules involved, the stochastic approach should be more cor- 
rect in describing the real system: we study the agreement between the two approaches by 
exploring the system parameter space. We address the problem by proposing a simplified 
version of the model that allows analytical treatment, and by performing numerical simu- 
lations for the full model. We observed optimal agreement between the stochastic and the 
deterministic description of the circuit in a large range of parameters, but some substantial 
differences arise in at least two cases: 1) when the deterministic system is in the proximity 
of a transition from a monostable to a bistable configuration, and 2) when bistability (in the 
deterministic system) is "masked" in the stochastic system by the distribution tails. The 
approach provides interesting estimates of the optimal number of molecules involved in the 
toggle. Our discussion of the points of strengths, potentiality and weakness of the chemical 
master equation in systems biology and the differences with respect to deterministic model- 
ing are leveraged in order to provide useful advice for both the bioinformatician practitioner 
and the theoretical scientist. 

Keywords: Toggle switch, chemical master equation, stochastic versus deterministic, bio- 
chemical reactions 

Introduction 

Complex cellular responses are often modeled as switching between phenotype states, and 
despite the large body of deterministic studies and the increasing work aimed to elucidate 
the effect of intrinsic and extrinsic noise in such systems, some aspects still remain unclear. 
Molecular noise, which arises from the randomness of the discrete events in the cell (for ex- 
ample DNA mutations and repair) and experimental studies have reported the presence of 
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stochastic mechanisms in cellular processes such as gene expression [50], [51], [35] . decisions 
of the cell fate [28J, and circadian oscillations [22] • Particularly, low copy numbers of im- 
portant cellular components and molecules give rise to stochasticity in gene expression and 
protein synthesis, and it is a fundamental aspect to be taken into account for studying such 
biochemical models [HE]- In this paper, we consider a simplified circuit that is known to 
govern a fundamental step during the eukaryotic cell cycle that defines cell fate, previously 
studied by means of a deterministic modeling approach [3]. Let set the scene by reminding 
that " all models are wrong, but some are useful" (said by George Edward Pelham Box, who 
was the son-in-law of Ronald Fisher). Biologists make use of qualitative models through 
graphs; quantitative modeling in biochemistry has been mainly based on the Law of Mass 
Action which has been used to frame the entire kinetic modeling of biochemical reactions 
for individual enzymes and for enzymatic reaction network systems [32]. The state of the 
system at any particular instant is therefore regarded as a vector (or list) of amounts or 
concentrations and the changes in amount or concentration are assumed to occur by a con- 
tinuous and deterministic process that is computed using the ordinary differential equation 
(ODE) approach. However, the theory based on the Law of Mass Action does not consider 
the effect of fluctuations. If the concentrations of the molecules is not great (on the order of 
Avogadro's number) we cannot ignore fluctuations. Moreover, biological systems also show 
heterogeneity which occurs as a phenotypic consequence for a cell population given stochas- 
tic single-cell dynamics (when the population is not isogenic and in the same conditions). 
From a practical point of view, for concentrations greater than about 10 nM, we are safe 
using ODEs; considering a cell with a volume of 10~ 13 liters this corresponds to thousands 
of molecules that, under poissonian hypothesis, has an uncertainty in the order of 1%. If 
the total number of molecules of any particular substance, say, a transcription factor, is 
less than 1,000, then a stochastic differential equation or a Monte Carlo model would be 
more appropriate. Similarly to the deterministic case, only simple systems are analytically 
tractable in the stochastic approach, i.e. the full probability distribution for the state of the 
biological system over time can be calculated explicitly, becoming computationally infeasi- 
ble for systems with distinct processes operating on different timescales. An active area of 
research is represented by development of approximate stochastic simulation algorithms. As 
commented recently by Wilkinson the difference between an approximate and exact model is 
usually remarkably less than the difference between the exact model and the real biological 
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process [7J . Given we can see this either as an unsatisfactorily state of art or as a promising 
advancement, we can summarise the methodological approaches as following. Biochemical 
networks have been modeled using differential equations when considering continuous vari- 
ables changing with time, or stochastic differential equations (SDE) for the trajectories of 
continuous random variables changing with time, and using the Gillespie algorithm that 
could be thought as an algorithm for the trajectory of discrete random variables changing 
with time. For a random variable changing with time one can either characterize it by its 
stochastic trajectories, as by the SDE and the Gillespie algorithm, or one can characterize 
its probability distribution as a function of time. The corresponding equation for the SDE 
is the Fokker-Planck equation, and the corresponding equation for the Gillespie algorithm 
is called the chemical master equation (CME) [33] . 

Therefore, the chemical master equation is a Markov chain in the limit where time is 
continuous and it could be thought as the mesoscopic version of the Law of Mass Action, 
i.e. it extends the Law of Mass Action to the mesoscopic chemistry and biochemistry, see 
for example [44, 45]. 

Here we compare the results of a stochastic versus deterministic analysis of a microRNA- 
protein toggle switch involved in tumorigenesis with the aim of identifying the most mean- 
ingful amount of information to discriminate cancer and healthy states. We show that the 
stochastic counterpart of such deterministic model has many commonalities with the deter- 
ministic one, but some differences arise, in particular regarding the number of stable states 
that can be explored by the system. The disagreement between the stochastic and determin- 
istic description is observed in a "ghost" effect caused by the proximity to a deterministic 
bifurcation [17], and in a somehow opposite situation, in which the variance of the stable 
point can mask the detection of the second peak in the stationary distribution. In this paper 
we perform a numerical study of the complete two-dimensional model, but we consider also 
a simplified, biologically meaningful, version of the model that allows to calculate an exact 
solution, with a numerical characterization of the parameter ranges in which the two systems 
produce qualitatively similar results. 
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I. PROPERTIES OF A MICRORNA TOGGLE SWITCH 



Oncogenes and tumor-suppressor genes are two pivotal factors in tumorigenesis. Recent 
evidences indicate that MicroRNAs (miRNAs) can function as tumor suppressors and onco- 
genes, and these miRNAs are referred to as oncomirs. MiRNAs are small, non-coding RNAs 
that modulate the expression of target mRNAs. The biogenesis pathway of miRNAs in 
animals was elucidated by [24]. MiRNAs undergo substantial processing since the nuclear 
transcription where two proteins play an essential role: Drosha and Dicer. Most of miRNA 
are first processed into pre- miRNA by Drosha. After exportated to the cytoplasm, the pre- 
miRNA is processed by Dicer into a small double strand RNA (dsRNA) called the miRNA: 
miRNA duplex. The active strand, which is the mature miRNA is incorporated into the 
RISC and binds to the target mRNA, whereas the inactive strand is ejected and degraded. 
In normal tissue, proper regulation of miRNAs maintains a normal rate of development, cell 
growth, proliferation, differentiation and apoptosis. Tumorigenesis can be observed when 
the target gene is an oncogene, and the loss of the miRNA, which functions as a tumor sup- 
pressor, might lead to a high expression level of the oncoprotein. When a miRNA functions 
as an oncogene, its constitutive amplification or overexpression could cause repression of its 
target gene, which has a role of tumor suppressor gene, thus, in this situation, cell is likely 
to enter tumorigenesis. MiRNAs are often part of toggle switches, with important examples 
are gene pairs built with oncogenes and tumour suppressor genes [5J E] • Here we focus on 
the amplification of 13q31-q32, which is the locus of the the miR-17-92. The miR-17-92 
cluster forms a bistable switch with Myc and the E2F proteins (3j [9J [10] . The oncogene Myc 
regulates an estimated 10% to 15% of genes in the human genome, while the disregulated 
function of Myc is one of the most common abnormalities in human malignancy [111 [13]. The 
other component of the toggle is the E2F family of transcription factors, including E2F1, 
E2F2 and E2F3, all driving the mammalian cell cycle progression from Gl into S phase. 
High levels of E2Fs, E2F1 in particular, can induce apoptosis in response to DNA damage. 
The toggle also interacts with dozens of genes (see figure [I] depicts a portion), particularly 
with the Rb and other key cell cycle players. A summary of the experiments perturbing 
miRNA/Myc/E2F and E2F/RB behaviours have suggested the following: 

• The Rb/E2F toggle switch is OFF when RB inhibits E2F, i.e. stopping cell prolif- 
eration; it is ON when E2F prevails and induces proliferation. Once turned ON by 
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sufficient stimulation, E2F can memorize and maintain this ON state independently 
of continuous serum stimulation. 



The proteins E2F and Myc facilitate the expression of each other and the E2F protein 
induces the expression of its own gene (positive feedback loop). They also induce the 
transcription of microRNA- 17-92 which in turn inhibits both E2F and Myc (negative 
feedback loop) 



Moreover, the increasing levels of E2F or Myc drive the sequence of cellular states, namely, 
quiescence, cell proliferation (cancer) or cell death (apoptosis). 



Myc 



E2F 



Rb 



<(miR-17-92)> 



Cdc25A* 



cycD/cdk4 



cycE/cdk2 



p27 



FIG. 1: The E2F-MYC-miR-17-92 toggle switch with its biochemical environment 



Although there is increasing amount of research on cell cycle regulation, the mathematical 
description of even a minimal portion of the E2F, Myc and miR-17-92 toggle switch is far 
from trivial. Aguda and collaborators [3] have developed a deterministic model, which 
reduces the full biochemical network of the toggle switch to a protein (representing the 
E2F-Myc compound) and the microRNA- 17-92 cluster (seen as a single element). 

It is a 2-dimensional open system, in which p represents the E2f-myc complex and m the 
miRNA cluster: thus no mass action law holds, i.e. p + m ^ N . The dynamics of p and m 
concentrations are described by eq. [T] and [2] 

ki ■ p 2 



p = a + 



T 1 + p 2 + T 2 • m 
m = (3 + k2 ■ p — 7 • m 



— 5 ■ p 



(1) 
(2) 
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The system can be rewritten in an adimensional form as follows: 

ed> = a + — Z — — — 6 (3) 

fi = 1 + (f) - fj, (4) 

Where the parameters are: a' = j^oi, k = ^l 1 , T[ = T' 2 = ip-r 2 , e = | and the change 

of variables is: = /i = and r = jt. 

In the original model [3] , the rate of protein synthesis is not a function of the instantaneous 
concentration (as assumed in eq|3] ) but rather of its concentration at some time A in the 
past: 

1 / k[<p(r-A)] 2 

e ^ a+ P 1 + [0(,-A)F + p 2 .Mr-A) - 0(r) ' (5) 

We will not consider such delay in our stochastic realization of the model, since it would 
increase system dimensionality and it does not seem necessary to obtain the features we 
want to characterize. 

The steady state can be studied in the nondimensionalized system and, therefore, the 
conditions on the parameters for the existence of multiple steady states. In the resulting 
cubic equation: 

the necessary (but not sufficient) conditions for the existence of 3 steady states (and thus a 
bistable system) are: 

(7) 



(la - *) < cl < (l + j|) 



II. THE STOCHASTIC MODELING APPROACH 

The system represented by equations [T] and [2] can studied as a stochastic system through 
the Chemical Master Equation (CME) approach [4j. The resulting CME has two variables, 
the number of p and m molecules, labeled as n and m. The temporal evolution in the 
probability, p n ,m(t), to have n and m molecules at time t is described by the following 
equation: 

p n , m = (E n - l)r n p nm + (E- 1 - l)g n p nm + (E m - l)r m p nm + (E m x - l)g m p nm (8) 
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This CME is derived under the condition of a one-step Poisson process, E and E -1 are the 
forward and backward step operators, g and r the generation and recombination terms for 
the n and m variables, as shown in superscripts. 

The two generation and recombination terms associated with the n and m variables are 
respectively: 

9 n = a+ r k \'f ; r n = 5-n (9) 

g m = /3 + k 2 -n; r m = 7 • m (10) 

We remark that the molecule influxes into the system (represented by the a and j3 terms) 
could be included in different ways in the stochastic equations, since in the deterministic 
equations they represent a sort of "mean field" value. As an example, molecules could be 
added in bursts with specific time distributions, that do not appear in the macroscopic 
continuous deterministic equations. We will consider the simplest approach, but the choice 
of different influx patterns should deserve further investigation. 

A. The one-dimensional model 

We can reduce the problem from two to one dimension, by considering a different time 
scale for the two reactions (in particular considering m^> p) and thus considering the steady 
state solution for the m: 

m = — p + k ■ p (11) 

7 

As a consequence we obtain a deterministic equation for the p only: 

Where V = and T" = r x + The stochastic equation for p n is thus as follows: 

p n = (E - l)r n ■ p n + (E- 1 - l)g n ■ p n (13) 

9n = a+ r + r».n + n* ; r ' n = 5 - n (14) 

A general solution can be obtained 

N ,. . N j fci-(i-l) 2 

n s _ TT gi! ~ X ) . „ _ IT r>+r»-(i-i)+(i-i)^ ( , 

1=1 v ' 1=1 
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PO = - ; ^ N z^v (16) 



with an adequate normalization factor imposed on p^: 

1 

i + EillllillPn! 

We remark that the system is open, thus in theory N is not fixed, but we can truncate 
the product to a sufficiently high value of N obtaining a good approximation of the whole 
distribution. This one-dimensional system (for which an analytical solution can be obtained) 
will be compared to numerical simulations of the exact one-dimensional and two-dimensional 
systems. 



III. MODEL ANALYSIS 



A. The stationary distribution 

The one-dimensional model can show monomodal as well as bimodal stationary distribu- 
tions, depending on the parameters considered. As an example, we obtain bistability with 
a set of parameters as shown in Fig. [2] 

0.014 

0.012 

0.010 

s; 0.008 
a 

g 0.006 

s 

0.004 
0.002 
0.000 

50 100 150 200 250 300 

U of molecules 

FIG. 2: The stationary distribution for the one-dimensional space, obtained using the following 
parameters: a = 1.68(molecule/h), (3 = 0.202(molecule/h), 5 = 0.2(/i -1 ), 7 = 0.2(/i _1 ), Ti = 
10300(molecule 2 ) , T2 = 1006 (mol ecul e) , k\ = 90(molecule/h) and ki = 0.05(/i~l). 

Thus the qualitative features of the two-dimensional deterministic model (i.e. the 
possibility of being bistable depending on the parameter range) are recovered for the one- 
dimensional approximation of the stochastic system. Also the two-dimensional stochastic 
system shows bistability for the same parameters, and they are in optimal agreement for a 
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probability distibution 




20 40 SO SO 100 120 140 160 

# of p molecules 

FIG. 3: Comparison between the deterministic solution (bottom) and the stationary distribution 
(top) for the parameter set as in Table [IJ case 6. 

range of parameters in which the rh ^> p condition holds 

We also observe some remarkable differences between the deterministic and the stochastic 
models: there are regions in parameter space in which the deterministic approach shows only 
one stable state, but in the stochastic system two maxima in the stationary distribution are 
observed (see Fig. [3J. This difference can be explained qualitatively as follows: for the 
deterministic system, there are parameter values for which the system is monostable but 
very close to the "transition point" in which the system becomes bistable. It is known that 
in these situations a "ghost" remains in the region where the stable point has disappeared 
|17j . for which the systems dynamics has a sensible slowing down (i.e. when the system 
is close to the disappeared fixed point, it remains "trapped" for a longer time close to it, 
in comparison with other regions). This behaviour results in the presence of a peak in the 
stationary distribution of the corresponding stochastic systems, that thus remains bistable 
also when the deterministic system is not anymore. 

Another difference is observed: for some parameter values the deterministic system is 
bistable, but the stochastic distribution shows a clear peak for the maximum with the largest 
basin of attraction and the smaller peak results "masked" by the tail of the distribution 
around the first peak (see Fig. [4]), thus resulting in a monomodal distribution with a long 
tail. In practice, the highest state behaves like a sort of metastable state, since the states 
of the system with a high protein level are visited only occasionally. 
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FIG. 4: Comparison between the deterministic solution (bottom) and the stationary distribution 
(top) for the parameter set as in Table [IJ case 4. 

B. Numerical analysis 

Here we implemented numerical methods to find the stationary distribution of a CME. 
The most accurate is the Kernel resolution method: given the complete transition matrix of 
the system, it is possible to solve numerically the eigenvalue problem, obtaining the correct 
stationary distribution. This method, in this case, has a serious drawback: the system is 
of non-finite size, preventing a complete enumeration of the possible states. Even with a 
truncation, the system size rises in a dramatic way: the state space for a bidimensional 
system is of order N 2 if N is the truncation limit, and thus the respective transition matrix 
is of order iV 4 . This means that even for a relatively small system (with a few hundred of 
molecules) the matrix size explodes well beyond the computational limits. The only feasible 
resolution strategy is a massive exploration of state space by Montecarlo methods, in which 
single trajectories of the system are simulated: performing this simulations long enough for 
several times allows to estimate the stationary distribution. 

The Montecarlo method we chose is a modified version of the SSA algorithm (also known 
as the Gillespie algorithm) named logarithmic direct method [201 IM]. which is a statistically 
correct simulation of an ergodic Markov system. It is not the fastest algorithm available, as 
compared to other methods like the next-reaction or the r-leap method, but it produces a 
correct estimation of the statistical dispersion of the final state. 

For each parameter set we performed 10 simulations for about 10 6 — 10 7 iteration steps 
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TABLE I: Table of the parameter sets for the cases considered. 



Par 


Case 1 


Case 2 


Case 3 


Case 4 


a (molecule/h) 


1.0 


1.68 


1.0 


20.0 


5 (h- 1 ) 


1.0 


0.20 


0.09 


1.19 


8 (molecule/h) 


1.0 


0.202 


0.0 


1.0 


1 V 1 


100.0 


0.20 


10.0 


1.0 


k\ (molecule/h) 


30.0 


90.0 


12.5 


230.0 


k 2 (h-1) 


100.0 


0.05 


10.0 


1.0 


Fi (molecule 2 ) 


60.0 


10300.0 


(72.8) 2 


(110.0) 2 


F2 (molecule) 


10.0 


1006.0 


10.0 


10.0 




U of molecules 




tt of molecules 



FIG. 5: Case of good agreement between the theoretical and obtained distribution (see Tab. |TJ 
case 1). Left: one-dimensional system, right: two-dimensional system. The thin black line is 



the theoretical distribution obtained from Eq. 15 The thick dark grey line is the average of the 
various simulations, while the grey and light grey areas represent the range of one and two standard 
deviations from the average distribution. 

each. The multiple simulations were averaged together for a better estimation of the sta- 
tionary distribution, and they allowed also an estimation of the variance over this average 
distribution. 

In the following we discuss four cases that describe the system behaviour for different 
parameter settings, shown in Table |TJ 
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FIG. 6: Case of poor agreement between the theoretical and obtained distribution (see Tab. 
case 2). Left: one-dimensional system, right: two-dimensional system. The thin black line is the 



theoretical distribution obtained from Eq. 15 The thick dark grey line is the average of the 
various simulation, while the grey and light grey areas represent the range of one and two standard 
deviations from the average distribution. 

In case 1, we have a system in which the hypothesis of a time-scale separation between 
m and p is strongly satisfied. The simulation was performed up to a time limit of 10 3 : we 
can see how the two resulting distributions are in good agreement with the theoretical one 
(see Fig. [5]), with the regions of higher variance of the histogram around the maxima and 
minima of the distribution. 

In case 2, the time-scale separation assumption does not hold, due to the very low value of 
7 and fo: even if this condition doesn't guarantee that the stationary state will be different 
from the approximate one-dimensional solution, with this set of parameters we can see a 
huge difference between the two distributions (Fig. [6]). 

In case 3, as defined before, we observe a "ghost" in which, even if a deterministic stable 
state does not exist, we can clearly see a second peak in the distribution (Fig. [7]). In this 
system the time-scale separation assumption holds, and we can see how both distributions 
show similar features. 

In this final case (Tab. |TJ case 4, Fig. [8J we can see another effect, in which the peak 
related to a deterministic stable state is masked by the tail of the stronger peak, becoming 
just a fat tail. Even without a strong time-scale separation for the m and p variables, we 
can see how both systems give a very similar response, evidencing that this effect is very 
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FIG. 7: Case 3, "ghost effect": only the biggest peak comes from a deterministic stable point. Left: 
one-dimensional system, right: two-dimensional system. The thick dark gray line is the average 
of the various simulation, while the gray and light gray areas represent the range of one and two 
standard deviations from the average distribution. 





ft of molecules 



tt of molecules 



FIG. 8: Case 4, peak masking effect (parameters as in Tab. | case 4). The deterministic system 
has two stable points, but only the peak related to the smallest stable point (with the largest basin 
of attraction) is visible. Left: one-dimensional system, right: two-dimensional system. 

robust. Increasing the 7 and &2 values does not affect the distribution as long as their ratio 
is kept constant. Note that while there are several computational tools for discrete-state 
Markov processes such as PRISM [SU], APNNtoolbox 02], SHARPE [H], or Mobius [H], 
there is very little for CMTC (see for instance [4"I]). Different modeling approaches for toggle 
switches do exists in the area of formal methods (see for example |2U 122]). 
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IV. DISCUSSION AND CONCLUSION 



We have studied a stochastic version of a biochemical circuit that is supposed to be 
involved in cell cycle control, with implications for the onset of severe diseases such as cancer, 
consisting of a gene cluster (Myc-E2F) and a miRNA cluster (mir-17-92). This cluster 
has been reported in very large number of cancer types: particularly in different types of 
lymphomas, glioma, non-small cell lung cancer, bladder cancer, squamous-cell carcinoma of 
the head and neck, peripheral nerve sheath tumor, malignant fibrous histiocytoma, alveolar 
rhabdomyosarcoma, liposarcoma and colon carcinomas. This huge variety of cancer stresses 
the centrality of this toggle switch and suggests that advancement in modeling this toggle 
could lead to insights into differences between these cancers. This aim is still far but we 
are delighted to report that our modeling approach shows important results inching to that 
direction. First of all, many features are recovered as observed for the deterministic version 
of the same system, also by means of a further approximation that reduces the system to 
an unique variable: in this case the system can be treated analytically, and compared to the 
one- and two-dimensional numerical simulations. 

The stochastic approach, that is the exact approach when the number of molecules in- 
volved is low, shows a different behaviour than the deterministic one in two situations we 
have observed. It is noteworthy that the number of molecules involved shows some agree- 
ment with the estimates by [2S] and by [23] for other miRNA-systems (see also [27]). The 
cell volume is assumed 10 -13 liters, then 1 nM =100 molecules. 

First, bistability in the stochastic system (namely, the possibility of having two stable 
states, one associated to a resting and the other to a proliferative cell state) is observed also 
in situations in which the corresponding deterministic system is monostable, and this can 
be explained by the presence of a "ghost" state in the deterministic system that is strong 
enough to produce a second peak in the stationary distribution of the stochastic model. 

Secondly, there are situations in which the peak for the stochastic distribution related to 
the highest level of expression (with parameter values for which the deterministic system is 
bistable) is masked by the tail of the distribution of the lowest-expression maximum (that 
is related to the largest basin of attraction in the deterministic model), making the "prolif- 
erative state" appear almost as a scarcely visited metastable state. This is an interesting 
behaviour, that should be further investigated in real experimental data of protein concentra- 
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tion and gene expression related to the biochemical circuit considered. The "metastable" and 
the "fully" bimodal distributions could be associated to healthy and tumoral cell states re- 
spectively, because the highest "proliferative" state has different properties in the two cases. 
From a biological point of view such state, being associated to a dysregulated, disease-related 
conditions, could actually represent a compendium of several dysregulated states. 

We argue that the deterministic approach to this biochemical circuit is not capable to 
characterize it completely, and the stochastic approach appears more informative: further 
features unique to the stochastic model could be obtained by considering different time 
patterns for the molecular influxes to the system, and this point in our opinion should 
deserve more investigation in a future work. MicroRNAs (miRNAs) express differently in 
normal and cancerous tissues and thus are regarded as potent cancer biomarkers for early 
diagnosis. We believe that the potential use of oncomirs in cancer diagnosis, therapies and 
prognosis will benefit accurate cancer mathematical models. 

Given that MiR-17-5p seems to act as both oncogene and tumor suppressor through de- 
creasing the expression levels of anti-proliferative genes and proliferative genes, this behavior 
is suggestive of a cell type dependent toggle switch. Therefore fitting of experimental data 
could provide insights into differences among cancer types and on which cell type is behaving 
differently. 
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